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Abstract 

Using Monte Carlo simulations we investigate some new aspects of the phase diagram and the 
behavior of the diffusion coefficient in an associating lattice gas (ALG) model on different regions 
of the phase diagram. The ALG model combines a two dimensional lattice gas where particles 
interact through a soft core potential and orientational degrees of freedom. The competition 
between soft core potential and directional attractive forces results in a high density liquid phase, 
a low density liquid phase, and a gas phase. Besides anomalies in the behavior of the density 
with the temperature at constant pressure and of the diffusion coefficient with density at constant 
temperature are also found. The two liquid phases are separated by a coexistence line that ends in 
a bicritical point. The low density liquid phase is separated from the gas phase by a coexistence 
line that ends in tricritical point. The bicritical and tricritical points are linked by a critical A-line. 
The high density liquid phase and the fluid phases are separated by a second r critical line. We 
then investigate how the diffusion coefficient behaves on different regions of the chemical potential- 
temperature phase diagram. We find that diffusivity undergoes two types of dynamic transitions: 
a fragile-to-strong transition when the critical A-line is crossed by decreasing the temperature at 
a constant chemical potential; and a strong-to-strong transition when the r-critical line is crossed 
by decreasing the temperature at a constant chemical potential. 

PACS numbers: 64.70.Pf, 82.70.Dd, 83.10.Rs, 61.20.Ja 
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I. INTRODUCTION 



The study of the properties of supercooled water is motivated by its well known anomalous 
thermodynamic behavior. Besides the density anomaly, the response functions for water 
appear to diverge at a singular temperature = 228 K [1]. This apparent divergence of 
the response functions led to the hypotheses of the existence of liquid polimorphism and 
of a second critical point, at = 228 K [2]. In spite of the enormous attention given to 
this possible singularity, as well as to the many other anomalies, no unique explanation has 
yet been established. The hypothetical singular point is hidden below the homogeneous 
nucleation temperature Th = 235 in an experimentally inaccessible temperature range 
for bulk supercooled water. This rules out direct experimental investigation of this region 
in order to confirm the existence of liquid-liquid coexistence. In order to circumvent this 
difficulty, it has been proposed, recently, that a dynamic crossover of the transport properties 
such as the self-diffusion constant, D, and the viscosity, rj, at temperatures above T^, would 
indicate the presence crossover has also been associated 

with liquid-liquid transitions in silicon 6| and in non-tetrahedral liquids Q]- 

The basic surmise behind the link between the dynamic crossover and the presence of a 
second critical point goes as follows. The liquid-liquid coexistence line that separates two 
liquid phases terminates at a critical point. Beyond this point, at which the response func- 
tions diverge, one finds lines of maxima of these functions which assymptotically approach 
the critical point. This extension of the first-order phase boundary into the one-phase region 
is the Widom line at T^(P). Even though this line does not exhibit any thermodynamic 
transition, experiments on water show that the specific heat, shear viscosity and thermal 



diffusivity 



exhibit a peak when crossing the Widom line. In particular, Maruyama et. 



al conducted experiments in nanopores (to avoid homogeneous nucleation) at ambient 
pressure that present a peak at the constant pressure specific heat at Tcp = 227 K. This 
temperature coincides (within the experimental error bar) with that one temperature ob- 
tained by Xu et. al 1^, Tcp = 225 K for the location of the dynamic crossover suggesting 
that this crossover occurs at the Widom line, confirming the presence of the second critical 
point. 

Unfortunately, the presence of a peak in the specific heat in a certain region of the 
pressure-temperature phase-diagram is not exclusivity of Widom lines. For instance, in 
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glassformers an abrupt heat capacity drop is observed when ergodicity is broken. This 
change can happen very sharply in the case of fragile liquids or it may take tens of degrees in 



the case of strong liquids. Examples of fragile liquids are to 



covalent and network forming systems are strong liquids 



uene and metallic systems, while 
In the last case, the increase in 
the specific heat can be simply a smeared peak, located above the melting temperature, T^, 
like in the case of Si02 and of BeF2 Q Q Q Q Q H- 

In the case of strong liquids, it 
is also possible to observe a weak transition at a temperature between the glass transition 
temperature, Tg, and the melting temperature, T^. This peak in the specific heat curve 
occurs in the tail of a A thermodynamic transition where there is a little heat capacity to 
loose [3| • Example of such strong liquids are the tetrahedral bonded liquids such as water. 
Si and Ge. This implies that observing a fragile-to-strong crossover in a region where the 
specific heat grows does not univocally imply the presence of a critical point. An interesting 
question, however, would be: does the presence of criticality result in a fragile-to-strong 
crossover? 

In order to address this point, in this paper we analyze a model that exhibits two different 
critical lines and we explore what happens with the dynamics close to these line, in order to 
test if a fragil-to-strong transition would be a signature for criticality. The present model is 
an Associating Lattice Gas (Henriques and Barbosa) that corresponds to a lattice gas with 
hydrogen bonds representede through ice variables. A competition between the filling up 
of the lattice and the formation of an open four-bonded orientational structure is naturally 
introduced in terms of the ice bonding variables, and no ad hoc addition of density or bond 
strength variations is needed. Besides the gas phase and as a result of this competition, the 
model exhibits two liquid phases that bare resemblance to the two liquid phases predicted 
for water, corresponding to a low density liquid phase and a high density liquid phase. 
Moreover, it has both the diffusion and the density anomalies present in water 19 1. 

Here, the model phase diagram is reviewed and analyzed for the presence of dynamic 
transitions. Two new critical lines were found beyond the liquid-liquid coexsitence line. We 
searched for fragile to strong transitions in the proximity of these two lines. Comparison 
between the behaviors of the specific heat and of the diffusion constant in these regions 
may help in understanding if the type of dynamic transition observed in confined water 
necessarily means the presence of criticality. 

The remaining of this article goes as follows. In sec. II, the lattice model is reviewed, for 
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clarity. In sec. Ill, results for the chemical potential-temperature phase-diagram are shown 
and discussed. Our investigation of diffusion is presented in sec. IV. Sec V is a final section 
of conclusions. 



II. THE MODEL 



We consider a two-dimensional lattice gas model of size L on a triangular lattice as 



introduced by Henriques and Barbosa 20|]. In this model, particles are represented by an 
occupational variable, 0"^, which assumes the value cr^ = 0, if the site is empty, or cr^ = 1, if 
the site is full, and six orientational variables, rf-, that represent the different orientations 
that the particle might exhibit. If two neighboring sites have complementary orientations, 
a hydrogen bond is formed. Four bonding variables are the ice bonding arms: two donors, 
with rf- = 1, and two acceptors, with rf- = —1. The other two arms, with r;^ = 0, do not 
form bonds, and are taken always opposite to each other, as illustrated in Fig.([T]). There is 
no restriction for donor /acceptor arms positions, thus there are eighteen possible states for 
each occupied site. 

The Hamiltonian includes two contributions: an isotropic, van der Waals like interac- 
tion, while the second interaction depends on the orientational degrees of freedom. Two 
neighboring sites, i and k, with pointing arms A and B, form a hydrogen bond if the prod- 
uct between their orientational variables is given by rf-rl^ = —1, yielding an energy per site 
e = E/L"^ = -V. For a non-bonding pair of occupied sites, the energy per site is e = —v + 2u, 
for M > 0. In spite of the fact that each molecule may have six neighbors, only four hydrogen 
bonds per particle are allowed. The overall energy of the system is given by 



{i,k) {i,k) 

where (7^ = 0, 1 are occupation variables, rf- = 0, ±1 represent the arm state variables, the 
summation {i, k) is over neighboring sites. 

Comparing the energies of the model at zero temperature two liquid phases, a low density 
(LDL) and high density (HDL) liquid phase are found, besides the gas phase. Fig. ([2]) 
illustrates the HDL and LDL phases. For high values of the chemical potential the lattice is 
fully occupied ( density p = 1 ) and the energy per site is e = —3v + 2u. At lower values of 
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FIG. 1: Particles in the model: An occupied central site i and its six bond variables, rf, with 
A = 1, ..,6. If = no bond is formed in spite of the configuration of the arm of the neighbor 
site. If r/^ = ±1 and the neighbor's arm is = ^1, a bond is formed. Dashed lines represent a 
non bonding configuration, while the solid line represents a bonding configuration. 

the chemical potential, /i, the soft core repulsion becomes dominant, and the lattice becomes 
3/4 filled, with density p = 0.75 and energy per site e = —^v. Like every other lattice gas, 
the model exhibits a gas phase, at very low chemical potentials. 

At zero temperature, the grand potential per site, = $/L^, is given by 

<f){T ^0) ^ {n- f^Y^a,) ^ E - fiN . (2) 

i 

By equating the grand potential of different phases, we find that the high density phase 
(HDL) coexists with the low density phase (LDL) at the reduced chemical potential /7 = 
yu/f = —6+8u/v. The coexistence between the LDL and the gas phases occurs at /i/v = —2. 

The properties of the system at finite temperatures were obtained from Monte Carlo 
simulations in the grand canonical ensemble, through the Metropolis algorithm. We present 
a detailed study of the model system, for L — 30. Some finite size scaling analysis was 
also undertaken, when necessary. Interaction parameters were fixed at u/v — 1, which 
corresponds to "repulsive" van der Waals interaction. Reduced parameters are defined by 

P 

p = - 

V 

T = ^ (3) 

V 

= -. 

V 

(4) 




FIG. 2: In the high density phase, HDL, the lattice is fuh, and an energy punishment arises, 
because two inert arms point to filled sites. In the low density phase, LDL, the lattice is 3/4 filled 
and particles are distributed over the lattice in such a way that the inert arms point only to the 
empty sites. There is no energy punishment, in this case. 

Equilibrium transitions were investigated through analysis of the system specific heat. First- 
order transitions point were located from hysteresis. The constant volume specific heat 
was calculated from simulation data obtained at constant chemical potential through the 
expression 

1 



Cv 



keT^V 



m'),vT 



{H^P)IVT. 



liVT 



(5) 



adapted from 



2l| to the lattice, p is the density, V is the volume and 5X = X — {X) with 



III. THE PHASE DIAGRAM 



The chemical potential-temperature phase diagram of the model was partially analyzed 
in previous work 20], which focused on the coexistence lines between the low and high 



density liquids. In this paper the JI — T phase-diagram is complemented by the analysis of 
the region beyond the coexistence line. 

The complete Jl — T phase-diagram is illustrated in Fig. (JSj) and goes as follows. At low 
reduced chemical potentials, /I, for all reduced temperatures, T only the gas phase is present. 
As the reduced chemical potential increases a low density liquid phase appears. This phase 
coexists with the gas phase along a first-order transition line at /I = 'pgas-LOLiT). For 
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even higher reduced chemical potentials a high density liquid phase emerges. This phase 
coexists with the low density liquid phase at the first-order line n = IJ-ldl-hdl 

(T). But 

what happens at the end of the two first-order lines? 

In order to answer the question, we have examined the specific heat at constant volume, 
cy , as a function of temperature, for fixed values of /I, in two regions of the JI — T phase di- 
agram: between the two coexistence lines and above the LDL-HDL coexistence line. Fig.Q 
illustrates our results. For /I = 0, between the two coexistence lines, cy has a peak at a 
reduced temperature T = Ta ~ 0.79, suggesting the presence of criticality. Similar behavior 
was observed for every investigated chemical potential between the two coexistence lines, 
indicating the presence of a critical line. We called this line A and represented it in Fig. ([3]) 
through a dotted line and square symbols. Above the liquid-liquid coexistence line, for 
JI = 2.5, the specific heat, cy, displays also a peak at T = ^ 0.71. We have examined 
a range of chemical potentials above the LDL-HDL coexistence line. A line of maxima of 
these peaks, named r was added to the Jt — T phase diagram, as shown in in Fig. (dashed 
line and circles). 

In order to check the nature of the two lines, A and r, in Fig. (I3l), the specific heat Cy 
within these regions was computed for different lattice sizes (L=10, 20, 30, 40, 50, 80 and 
100). Fig. ([5]) illustrates the behavior of cy for /I = for various lattice sizes, showing a 
diverging peak, as L ^ oo. Fig. ([6]) shows cy for Ji = 2.5, for various lattice sizes. In this 
case, however, the peak increases mildly with L. 

The criticality of A and r was investigated by calculating the energy cumulant given by 

Fig-O illustrates the energy cumulant for /i = 0, showing the signature for criticality. Fig.([H]) 
illustrates the energy cumulant for fi = 2.5 that also indicates the presence of criticality. 

In the attempt to understand the differences between the two transitions, it is important 
to stablish what is the structural difference between the LDL, HDL, the high and the low 
densities fluid phases. To answer this question it is necessary to establish a measure of how 
structured is the liquid. We adopt the number of hydrogen bonds per particle, phb, and its 
correlation with particle density p as such a measure. Fig.Q shows that, as temperature is 
decreased towards the speciflc heat peak position, hydrogen bond density increases, while 
particle density decreases. This is indicative that bonding is accompanied by particles 
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FIG. 3: Phase diagram showing reduced chemical vs. reduced temperature. The diamonds repre- 
sent the Gas-LDL coexistence hne. The triangles indicate the LDL-HDL coexistence line. Td is 
the tricritical point Gas-LDL and Tc2 is the tricritical point LDL-HDL. The squares and circles are 
lines, obtained by the maximum in specific heat, that separates fluid phase from LDL and HDL 
phases, respectively. The zero temperature points, at /I = —2 and JI = 2, are exact. 

abandoning the lattice. On the other hand, in the case of the r line, both density and number 
of bonding particles increase as the temperature decrease towards the peak temperature, as 
shown in Fig. ffTOl) . Thus, in this case, bonding and lattice filling occur simultaneously. 
Fig.l fTTl) compares the behavior of bond density phb over the two transitions, while Figs. 
f[T^ and ([T3D illustrate the difference in density behavior. In the case of the A transition, 
as temperature decreases towards the LDL phase, the density, shown in Fig. f|T2l) . decreases 
drastically . At the transition, the system orders itself by forming hydrogen bonds and by 
releasing nonbonded particles. For temperatures below the transition, the density increases. 
As for the r line, density p, shown in Fig. (fT3|) .increases smoothly as temperature is lowered. 

A closer look of this idea is possible if one examines the behavior of the two densities on 
different sublattices. Inspection of Fig. ([2]) suggests dividing the lattice into four sublattices, 
as illustrated in Fig. lIT^ . Figs. lITS!) and (ITSl) display sublattice density variations and how 
the number of hydrogen bonds changes with temperature, in the critical region, for /I = 0. A 



9 




T 

FIG. 4: Specific heat at constant volume for Ji = and JI = 2.5. 

clear critical transition is seen, in which one sublattice becomes empty, while the other three 
get ordered. Figs.( fT7l) and ( fTSl) illustrate density variations and the number of hydrogen 
bonds of the four sublattices with temperature, in the r transition region, Ji = 2.5. In this 
case, both densities of the four sublattices change smoothly. 

Analysis of the sublattice data shows that A is an order-disorder transition in which, 
as the temperature is decreased from the disordered fluid phase, bonds are formed, while 
nonbonded sites become empty. This critical A line joins the two coexistence lines at a 
tricritical point, T^i, and a bicritical point, Tc2- 



IV. DYNAMICS 



In order to quantify mobility in supercooled liquids, the concept of fragility was introduced 



by Angell 22|. Analyzing relaxation as a function of temperature, liquids are classified as 
strong, when relaxation follows an Arrhenius law, or fragile, when the relaxation follows a 
non-Arrhenius law. Strong liquids present structure that is preserved when temperature is 
increased, whereas in fragile liquids this structure is easily broken, as temperature increases. 
Within the framework of the Adam-Gibbs theory Id], viscous liquids are described as 

10 
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FIG. 5: Specific heat at constant volume for different lattice sizes in the fluid/LDL region (// = 0). 
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FIG. 6: Specific heat at constant volume for different lattice sizes in the fluid/HDL region {fi = 2.5). 



being made of clusters that rearrange cooperatively in order to pass through the free energy 
barrier. Consequently, diffusion depends on this cooperative rearrangement of the clusters 
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FIG. 7: Energy cumulant for the fluid/LDL region (/x = 0). 
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FIG. 8: Energy cumulant for the fluid/HDL region (/x = 2.5). 



through equation 



(7) 



for the diffusion constant D. Here Dq and C are constants, A/x is the free energy barrier 
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FIG. 9: (Color on line) Density, specific heat and number of bonds per site versus reduced tem- 
perature for /I = 
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FIG. 11: Density of hydrogen bonds as a function of the reduced temperature. For Ji = 2.5 a 
smooth change in phi, is observed while for JI = has a more abrupt change characterizing a 
continuous transition. 
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FIG. 12: Density as a function of the reduced temperature, for /x = 0. 



which the clusters have to overcome. Sc is the configurational entropy, given by 



Tk 



T 



dT 



(8) 



that describes how the structure of the hquid changes with temperature. In Eq.(IH]) Tk 
is the Kauzmann temperature 22] (for which Sc{Tk) = 0) and Acp is the difference in 
specific heat between the crystal and the liquid configurations, at temperature T. If Sc, is 
temperature independent, the diffusion follows an Arrhenius law, the liquid is very structured 
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FIG. 13: Density as a function of the reduced temperature, for /I = 2.5. 
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FIG. 14: The lattice is divided into four sublattices 

and the system is a strong liquid. If the configurational entropy depends on temperature, 
Sc = ACpln T/Tk, Eq.(j7]) becomes a Vogel-Fulcher equation, the liquid is not structured 
and is classified as a fragile liquid. 

Now we investigate the dynamic properties on crossing the A and the r lines, at constant 
chemical potential (see Fig. ([3])), by analysing behavior of model diffusivity. In order to 
compute diffusion coefficient we first equilibrate the system at fixed chemical potential and 
temperature. In equilibrium this system has n particles. Starting from this equilibrium 
configuration at a time t = 0, each one of these n particles is allowed to move to an empty 
neighbor site randomly chosen. The movement is accepted if the total energy of the system 
is reduced, otherwise it is accepted with a probability exp{AE/kBT) where AE is the 
difference between the energy of the system after and before the movement. After repeating 
this procedure nt times, the mean square displacement per particle at a time t is given by 

(Ar(t)2) = ((r(t)-r(0)f) , (9) 

where r(0) is the particle position at the initial time and r{t) is the particle position at a time 
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FIG. 15: Density of each one of the four sublattices as a function of temperature for /i = 0. 
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FIG. 16: Hydrogen bonds per site of each one of the four sublattices as a function of temperature 
for Jl = 0. 



t. In Eq. ( TTTI) . the average is taken over all particles and over different initial configurations 
The diffusion coefficient is then obtained from Einstein's relation 

{Arm 



D = lim 



(10) 



t^co At 

Since the time is measured in Monte Carlo time steps and the distance in number of lattice 



distance, a dimensionless diffusion coefficient is defined as 



D = lim 



{Arm 



t^oo At 
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(11) 
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FIG. 17: Density of each one of the four sublattices as a function of temperature for /x = 2.5. 
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FIG. 18: Hydrogen bonds per site of each one of the four sublattices as a function of temperature 
for JI = 2.5. 

where r = r/a and a is the distance between two neighbor sites and t = t/tMc is the time 
in Monte Carlo steps. 



Figs. (19(a)) - (19(d)) illustrate the behavior of the diffusion constant D with the inverse 



of the reduced temperature 1/T, for fixed values of the reduced chemical potentials (/x = 
—0.25, 0, 0.5, 1) At higher temperatures, diffusivity follows a non-Arrhenius trend, namely 

y = A0 + A1X + A2X^ + A^x^ (12) 

indicating that the low density disordered fluid phase is a fragile liquid. At lower tempera- 
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(c)Diffusion in function of temperature for JI = 0.5. (d)Diffusion in function of temperature for 11=1. 

FIG. 19: (Color on line). The circles are diffusion coefficient measured in simulation, solid line is 
a cubic fit and dashed line is an exponential fit. In high temperatures system behaves as a fragile 
liquid following a non-Arrhenius law, while for low temperatures the system behaves like a strong 
liquid following an Arrhenius law. 



tures, diffusivity displays Arrhenius behavior, given by 



y = Aoexp[ 

X 



(13) 



thus characterizing the low density ordered hquid phase as a strong liquid. Ai are fitting 
parameters in both equations. 

This change in dynamics over the critical A line occurs because the liquid is structurally 
different on both sides of the critical line. In the low density disordered fluid phase, intersti- 
tial particles weaken the hydrogen bonds and disrupt the network, so particles can rearrange 
fast and the process of diffusion is not energy activated. In the LDL phase, the network is 
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(a)DifFusion in function of temperature for JI = 1.75. (b)DifFusion in function of temperature for /I = 1.85. 
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(c) Diffusion in function of temperature for /I = 1.95. 

FIG. 20: (Color on line). Circles are diffusion coefficient measured in simulation and solid and 
dashed lines are two different exponential fits. In the region of r line system behaves as a strong 
liquid in both sides. The difference between the two Arrhenius behaviors is the activation energy, 
that defines the slope of the curve. 



fully developed, resulting in an ordered liquid, in which particles are "trapped", increasing 
relaxation time and characterizing this phase as a strong liquid, in which an energy acti- 
vated diffusion process takes place. This is the dynamic transition observed when crossing 
a Widom line in ramp- like models QQ, 

which suggests that the dynamic transition is 
not linked with the type of line but with the structuring of the system if this happens with 
or without a thermodynamic phase-transition. The system becomes more organized, as can 
be seen from the drastic change in the density of the sublattices shown in Fig. ( fTSl) with a 
smooth change in the total density. 

Since the HDL is also a structured phase, in principle a fragile-strong transition in the 
dynamics of diffusion could also be expected on crossing the r line. However, this is not the 
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case. Figs. (20(a)) - (20(c)) illustrate the behavior of the diffusion constant as function of 
inverse temperature, 1/T, for fixed chemical potentials fl = 1.75,1.85 and 1.95. At higher 
temperatures and high chemical potentials (or equivalently high densities), the fluid phase 
has an Arrhenius behavior and so it is a strong liquid. At lower temperatures, the HDL 
phase also displays an Arrhenius behavior, and therefore is also a strong liquid. The HDL 
phase and high density fluid phases are both strong liquids that differ in the activation 
energy. In resume, when the system crosses the r line, we have a dynamic transition, and a 
strong-strong crossover is observed. In this case, the activation energy of the HDL phase is 
higher than the activation energy of the high density fluid phase, indicating that the HDL 
phase is more ordered than the high density fluid phase. Diffusion is lowered in the HDL 
phase because particles spend more time trying to rearrange, in comparison with the high 
density fluid phase. 

How can we explain the existence of a fragile-to-strong crossover on the critical A-line 
and a strong-to-strong transition on the r line? The answer is given by the structure of the 
liquid, described in the previous section. On crossing the A-line, the hydrogen-bonded net 
breaks down abruptly (see Fig. [T6l) . while the r-line is accompanied by a much smoother 
melting of the h-bond network (see Fig {18]). The difference in the change of structure on 
crossing one and the other lines is even clearer if one looks at sublattice densities. See 
Fig.dUD and Fig.([E 



V. CONCLUSIONS 



In this paper we have analyzed equilibrium and dynamic properties of the Associating 
Lattice Gas Model, a lattice gas with hydrogen bonds represented through ice variables. 
Competition between the filling up of the lattice and the formation of an open four-bonded 
orientational structure leads to the presence of two liquid phases and a gas phase. The 
coexistence lines between the LDL and the gas phases, and between the LDL and HDL 
phases are connected by a critical A-line. Besides the A-line, a second one, the r-line, also 
emerges from the LDL-HDL coexistence line. This line is also identified by a peak in the 
specific heat. 

The system undergoes two kinds of dynamic transitions: one fragile-to-strong (crossing 
the A-line) and a strong-to-strong (crossing the r-line). Both dynamic transitions are related 
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with structure of the system. In the fragile-to-strong case, the system change its structure 
drastically when cross A-line, from a disordered structure in high temperatures to an ordered 
one at low temperatures. In the r-line region, the change in structure is more subtle, but 

again the system undergoes a structural change between two ordered phases. 

Our results point out in the direction that criticality does not necessarily means fragil- 
strong transtion. This change is in fact related to the change of structure that in the present 
case appears in two very different forms. 
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